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ABSTRACT. We describe an approach for incorporating radiative transfer into 3D hydrodynamic 
cosmological simulations. The method, while approximate, allows for a self-consistent treatment of 
sclf-shiclding and shadowing, diffuse and point sources of radiation, and frequency dependent transfer. 
Applications include photodissociation, photoheating, and photoionization of the IGM. 



1. Introduction 



According to current thinking, the epoch of hydrogen reionization begins with the for- 
mation of the first massive stars in subgalactic objects at high redshifts z ~ 15 (e.g., 
Couchman & Rees 1986; Gnedin & Ostriker 1997; Haiman & Loeb 1997) and is es- 
sentially complete by z ~ 5 as required by the Gunn-Peterson test. Due to its higher 
ionization threshold and lower photoionization rate, helium reionization may be delayed 
to z ~ 3 or less (Meiksin & Madau 1993; Reimers et al. 1997). The essential reason for 
this delay is that it is not until z ~ 4 that the spectrum of the UV background has 
£h ' hardened sufficiently due the rising quasar population. 

Reionization has been studied in a spatially averaged fashion by many authors (e.g, 
Shapiro & Giroux 1987; Miralda-Escude & Ostriker 1990, 1992; Meiksin & Madau 1993; 
Giroux & Shapiro 1996; Haiman (these proceedings, and references therein.) These stud- 
ies reinforce Shapiro & Giroux's conclusion that quasars are insufficient to reionize the 
universe at z > 3.5 as required by the Gunn-Peterson measurement. (For a discussion 
of the uncertainties involved, see Meiksin & Madau 1993.) Barring exotica, such as 
decaying neutrinos, that leaves only stellar sources of reionization at higher redshifts. 
These stars will inevitably form in highly localized regions of space (i.e., galaxies.) Thus, 
cosmic reionization is intrinsically highly inhomogeneous. 

In this short contribution, we outline our intended approach to numerically simu- 
late inhomogeneous reionization in three spatial dimensions. We seek efficient numerical 
methods for evolving the ionizing radiation field in 3D which we may couple to hydrody- 
namic models of structure formation (e.g., Abel, Bryan & Norman, these proceedings.) 
We warn the reader that at this time, we have not proven that our method will be 
sufficiently accurate and robust to meet our goal, nor do we have any measurements of 
the computational cost. 



2. Description of the Problem 

Immediately after the first photoionizing sources turn on, one has a collection of isolated 
ionization zones (essentially HII regions) growing in an inhomogeneous, expanding uni- 
verse. Their number and evolution will depend on the source population, ionizing fluxes, 
and ambient conditions. Until the HII regions begin to interact (percolation), they can 
be solved as isolated cases. The homogeneous case has been solved by Shapiro (1986) for 
the case of quasar reionization. With no inhomogeneities to spoil spherical symmetry 
or emit appreciable ionizing recombination radiation, one can reduce the problem to 
an ODE for the radius of the I- front versus redshift. One avoids solving the radiative 
transfer equation by simply attenuating the ionizing flux by the 1/r 2 geometrical fac- 
tor and volumetric losses due to recombination to atomic levels n > 1. The peculiar 
velocity of the I- front is set by balancing, in the rest frame of the I- front, outgoing ion- 
izing photons and incoming neutrals. With these approximations, Shapiro found that 
for typical quasar luminosities, the I-front always expands supersonically with respect 
to both neutral and ionized media. Such I-fronts are known as weak R-type I-fronts (cf. 
Spitzer 1978). Thus hydrodynamic motions are unimportant, justifying their neglect a 
posteriori. Secondly, he found that the I-fronts never reach their equilibrium Stromgrcn 
radii in a Hubble time, although they do overlap completely by z = (but not by z = 5) 
for the observed number density of quasars. This is simply because the recombination 
time, which enters in the definition of the Stromgren radius, is longer than the Hubble 
time, and consequently the Stromgrcn radius is very large. 

The evolution of an HII region in an inhomogeneous medium (expanding or other- 
wise) in 3D is an unsolved problem — one that requires more powerful numerical meth- 
ods to solve in the general case. The principal complication is that localized density 
enhancements will both absorb the primary ionizing flux (shadowing), and isotropically 
emit ionizing radiation via recombinations to the ground state (diffuse radiation.) As a 
consequence, the ionizing radiation field becomes inhomogeneous and anisotropic. The 
presence of density enhancements (clouds) retard the expansion of the I-front relative 
to the homogenous case due to enhanced down-conversion of ionizing photons into non- 
ionizing Balmcr continuum photons. When the areal covering factor of opaque clouds 
approaches one, the I-front will become starved of ionizing photons and stop expanding. 
Well before this limit is reached, hydrodynamics will become important for two reasons. 
First, dense clouds inside the HII region will be photo-evaporated, leading to peculiar 
velocities of order the thermal sound speed ~ 30 km/ s (see Shapiro, these proceedings.) 
Secondly, when the I-front expansion rate falls below the sound speed in the HII region, 
it makes a transition to a weak D-type front (cf. Yorke 1986). When this occurs, the 
pressure difference between the ionized gas and the neutral gas drives a shock wave into 
the ambient medium. The subsequent expansion of the HII region is due to a combina- 
tion of hydrodynamic and radiative effects. A shell of material is accumulated between 
the inner I-front and the outer shock front, which separate from one another. 

The complexities of the percolation phase depend on whether hydrodynamic ef- 
fects are imporant or not. That in turn depends on the presence of density inhomo- 
geneities and whether the isolated HII regions become weak D-type before or after 
overlap. Shapiro & Giroux (1987) modeled the effects of density inhomogeneities on 



quasar-driven I-fronts phenomcnologically by adding a clumping factor to the Shapiro 

(1986) formalism. Physically, these clumps correspond to the Lya forest and Lyman 
limit systems. Meiksin & Madau (1993) calculate that the clouds may absorb as much 
as 50% of the UV radiation from point sources. Consequently, the clouds are a sub- 
stantial diffuse source of ionizing radiation (Haardt & Madau 1996.) Shapiro & Giroux 

(1987) showed that although the enhanced recombinations increase the UV background 
deficit relative to observed QSO source population, I-fronts remain supersonic (globally) 
until full overlap is achieved. 

However, if quasars fade, or if reionization is caused by far more numerous, less 
luminous stellar UV sources, hydrodynamics will become important as I-fronts become 
weak D-type as they approach their Stromgrcn radii. We desire numerical methods 
which can treat both kinds of circumstances. 



3. Basic Approach 

The basic elements of our method can be described very simply. In the next section 
we provide equations. We decompose the radiation field into point source and diffuse 
components: I v = If, ts + I^ l ff. The point source radiation field is attenuated along 
radial rays from each point source in the volume. Every cell is crossed by at least one 
ray. The number of photoionizations in frequency interval v, v + dv in each ray segment 
inside a cell is simply related to the decrease in If, ts along that segment. The total 
number of photoionizations in the cell is the sum over ray segments within the cell. The 
diffuse radiation field is computed by solving the angle-integrated moment equations 
with appropriate Eddington factor closure (Stone, Mihalas & Norman 1992). In the 
limit of small (compared to the horizon) volumes, the zeroth and first moment equations 
can be simplified and combined into a single nonlinear elliptic equation for J^ff = 
j- J I^ffdQ. This elliptic equation is discretized on a 3D cartesian mesh and may 
be solved using a variety of techniques, including multigrid relaxation and/or iterative 
sparse-banded matrix methods. 



4. Formalism 

4.1. Equation of Cosmological Radiative Transfer 

The equation of cosmological radiative transfer in comoving coordinates (cosmological, 
not fluid) is (Paschos, Mihalas & Norman 1998): 

\3h , n-V/„ H(t) 8I V 

— 5T + = (v- 3J„) = r\ v - xJ v (1) 

c at a c dv 

where I v = I(t, x, CI, v) is the monochromatic specific intensity of the radiation field, h 
is a unit vector along the direction of propagation of the ray; H(t) = a /a is the (time- 
dependent) Hubble constant, and a = is the ratio of cosmic scale factors between 
photon emission at frequency v and the present time t. The remaining variables have 
their traditional meanings (c.g, Mihalas 1978.) Equation (1) will be recognized as the 
standard equation of radiative transfer with two modifications: the denominator a in the 



second term, which accounts for the changes in path length along the ray due to cosmic 
expansion, and the third term, which accounts for cosmological redshift and dilution. 

In principle, one could solve equation (1) directly for the intensity at every point given 
r\ and x- However the high dimensionality of the problem (three positions, two angles, 
one frequency and time) not to mention the high spatial and angular resolution needed 
in cosmological simulations make this approach impractical for dynamic computations. 
Therefore we proceed through a sequence of well-motivated approximations which reduce 
the complexity to a tractable level. 

4.2. Local Approximation 

We begin by eliminating the cosmological terms and factors. That we can do this can be 
understood on simple physical grounds. Before the universe is reionized, it is opaque to 
H and He Lyman continuum photons. Consequently, ionizing sources are local to scales 
of interest, and not at cosmological distances. If our simulation box is of side length L 
and \ p is the photon mean free path, then by construction A p <I. The ratio of the 
third to the second terms in equation (1) is HLa/c -C 1, and hence the third term can 
safely be ignored. Now, let us consider the factor a in equation (1). For a photon which is 
emitted at time t on one side of the box and absorbed on the other side at time t + L/c, 
a = ( ~~f~ ) r * ~ l + r/L/ct = l + rjL/Ljj, where 77 is the logarithmic expansion rate of the 
universe (2/3 for il D = 1) and Lh is the Hubble horizon scale. For L <C Lh, a = 1, and 
v em = v. In practice, our dynamical timesteps are much longer than a photon crossing 
time. However, even in this case accuracy limits our dynamical timesteps such that 
Aa/a <C 1, and hence a = 1 in any given timestcp. Therefore, setting a = 1, equation 
(1) reduces to its standard, non-cosmological form: 

1 dl v 

--^r+n-WI v =r] v -XuI v (2) 
c at 

where now v is the instantaneous, comoving frequency. 

Now consider the case where X p » i, i.e., the simulation volume is optically thin to 
ionizing radiation as it would be after reionization. In that case ionizing sources are are 
either inside the box (local) or outside the box (nonlocal), or both. Local sources are 
treated as in the case above. Nonlocal sources sufficiently far from the box contribute 
to a nearly isotropic, homogeneous UV background (metagalactic flux). In this limit, 
equation (1) can be solved in an angle and spatially averaged fashion, and the cosmolog- 
ical term is not ignorable. This computation has been done by Haardt & Madau (1996) 
including emission from the observed quasar population, and absorption and re-emission 
by the Lya forest. The result of this calculation is J„(z) — the mean metagalactic in- 
tensity as a function of redshift. If the material in the simulation volume is optically 
thin everywhere, then the local radiation field is J„{z). Current simulations of the Lya 
forest (e.g., Zhang et al. 1997, 1998) employ this approximation. If, however, there exist 
opague regions within the simulation volume — for example, high column density Lyman 
alpha clouds — then J u no longer equals J* (z) locally, but must be computed using some 
approximation to equation (1) using J*(z) as a boundary condition. In this case, the 



cosmological terms are accounted for in the boundary conditions, and for box sizes much 
smaller than the horizon, ignorable within the box. 

4.3. Angular Moments 

While the radiation field due to local point sources is highly anisotropic, the diffuse 
radiation field should be more nearly isotropic since recombination radiation is emitted 
isotropically and absorbed by density enhancements which are well resolved on our 
computational grid. Thus we expect the angular structure in to be well described by 
its angular moments, the first three of which are defined as follows: J„ = j- § dVlI v \ K % v = 
§ dfln l I v ; K l J = ^ § dQn l n 3 I v . The radiation energy density, flux, and pressure 
tensor are related to these moments via the simple relations E u = ^-J u ,Fl = ^-H l u , 
and P l J = ^-K l J. 

Now, it is advantageous to work in a frame which is comoving with the fluid because 
in this frame the emission and absorption coefficients are isotropic (Mihalas & Mihalas 
1984; hereafter MM). Denoting comoving frame quantities with a subscript "o" , The first 
two moment equations of the radiation field are obtained by taking angular moments of 
equation (2), which yields temporarily suppressing the v subscript: 

) + V*^ + PpViVj = 4tt % - c Xo E , (3) 

where p is the fluid density and = ^ + v • V is the convective derivative. Equations 
(3) and (4) are MM equations (95.87) and (95.88) dropping acceleration terms. 



4.4. Quasi-Static Approximation 

Equations (3) and (4) can be simplified further when we realize we are interested in 
phenomena occuring on the fluid flow timescale, not the radiation flow timcscale. This 
quasi-static approximation amounts to throwing away terms which are always O(^) or 
higher. This must be done with care, and we refer the reader to MM for a thorough 
analysis. The result, for continuum radiation is: 

S7iFt = e v -ck v {t)E v (5) 

VjP n = - k -^Fl (6) 

where now it is understood all quantities are measure in the comoving (fluid) frame, 
e„ = A-nVv, and fc„ is the absorption coefficient (we have ignored scattering.) Equations 
(5) and (6) are called quasi-static because the only timescale which enters is through 
the material opacity k u (t), which evolves on a photoionization timescale. 



4.5. Closure Schemes 



Equations (5) and (6) are two equations in three unknowns: E, F and P. We will exper- 
iment with various approaches to closing the hierarchy of moment equations. The most 
general and accurate approach is the Eddington factor closure: = pjE v , where f l j 
is the tensor Eddington factor (e.g., Stone, Mihalas & Norman 1992). f % J contains all 
the angular information of the radiation field. If one knew /„(£!, x), say from solving 
equation (2), then can be computed via angular quadratures at every point x. How- 
ever, the goal is to avoid solving the angle-dependent transfer equation. One approach, 
which we call NEWS, is to compute l v only along rays parallel to grid lines, as well 
as along the principal diagonals. Since we use uniform cartesian meshes in cosmology, 
we simply evaluate the formal solution to eq. (2) on long characteristics. A second, less 
accurate approach, is to compute pi using geometric information about the location 
of principal emitters and absorbers. This approach would have to be calibrated against 
more exact methods. 

Yet more approximate, but perhaps adequate for our needs, is the diffusion approx- 
imation which states that F l v = —DV l E l/ , where D depends on the energy density and 
opacity through the relation: D = -r-X(E v ). The quantity A is called the flux limiter, 
and for optically thick media has a value of 1/3. When radiation propagates through op- 
tically thin media, it streams rather than diffuses. There are many functional forms for 
A which are taylored to the problem under study; we mention three which are prominent 
in the literature by Alme & Wilson (1974), Mincrbo (1978) and Levermore & Pomran- 
ing (1981). All depend on the quantity R — ^-7^7, which is used as a switch between 
the diffusion and free streaming limit. Whether any of these formulations will prove 
adequate for our application remains to be seen. 

Inserting the relation = —DV l E L/ into equation (5), we obtain 



where £s„ = — ^rS v , S: the source function. Since A = \{E U ,VE V ), equation (7) 
is a nonlinear, elliptic equation for E v , where the quantities and k v are functions 
of space and time. Solution requires the specification of boundary values on E v and its 
normal derivative. 

4.6. Frequency Reduction 

Finally, we consider reducing the frequency complexity For this, we employ the multi- 
group method (cf. MM, Ch. 6), in which the frequency spectrum is divided into a number 
of frequency groups. Since, in the first instance, we are are only interested in the pho- 
toionization of primordial gas, we need only consider three frequency groups above the 
ionization edges for HI, Hel and Hell at energies hv = 1, 1.809, and 4 Ryd, respectively. 
Defining the group average radiation energy density as: 



V i {—V i E v ) = k v {%„-E v ), 



(7) 




(8) 



the multigroup diffusion equations to be solved are: 

V l (DgV l Eg) =€g- C k„ E g , (9) 

J "s+i Du sjE u du _ J" ' 9+1 k v E v dv 

where D a = "?■„ and k„ = U L , the flux mean and absorption mean 

9 \ 9+1 VE„dv g 3+1 E v du ' ^ 

opacities, respectively. Because we don't know E v and VE V in advance, these mean opac- 
ities must be approximated. The approximations, while not unique, can be constructed 
to have various desirable properties, such as giving the exact energy or momentum 
absorbed within a group. We desire group means which conserve the total number of 
photoionizations within a group to a high degree of accuracy. Thus, we write: 

Vf J.. °^ v (~HI i _HeI i ,J?eJJ> \ " c ^g „ ,g; , .geJ , „HeII ^ „\ 
= / ~L v ^ +°V + a v > Z^T, < G v +(T " +(J v > 9 V^s+i -Vg) 

Jul hV ~i nv g+l/2 

(10) 

where vl is the Lyman limit, and af, is the photoionization cross section for species x. 
The angle average is defined formally as: 

<Q>g= / dvWg{v)Q{v)/{Vg +1 - Ug), (11) 

JVg 

where the weighting function W g {v) depends on the assumed spectral form (e.g., piece- 
wise constant) of E v within each frequency group. 

5. Coupling to Chemistry and Hydrodynamics 

Here we briefly discuss the solution strategy for the coupled radiation, matter, cos- 
mological fluid dynamical system. A method for efficiently computing nonequilibrium 
ionization, chemistry, and cooling coupled to mulitspecies cosmological hydrodynamics 
is described in Anninos et al. (1997). Applications of this method to first structure for- 
mation and the Lya forest are found in Abel et al. (1998) and Zhang et al. (1998), 
respectively. In the latter, photoionization of the IGM in the optically thin limit due to 
a metagalactic UV background is included. At the heart of the method is an implicit 
scheme for solving the kinetic equations for all the ionization states of H and He, as 
well as the reactants for Hi. Inside the main loop of the hydrodynamic computation, 
we subscycle on heating and cooling portion of the gas energy equation and the kinetic 
equations. The timestep during subcyclying is chosen such that the fractional abun- 
dances of species which dominate the cooling change by no more than 10%. Since all 
atomic and molecular cooling is proportional to the electron density, we find it sufficient 
to limit our chemistry timestep to At e = 0.1n e /h e . 

Radiative transfer merely changes when and where material gets photo- heated, ion- 
ized and dissociated. Thus, within the chemistry/cooling subcycle loop we also call the 
methods described above for computing the point source and diffuse radiation field. 
Since opacities change on the photoionization timcscalc which is reflected in the elec- 
tron fraction, we need not change our subcyle timestep criterion to include radiative 
transfer. 
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